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Abstract 

The starting point is the problem of finding the interaction energy 
of two coinciding homogeneous cubic charge distributions. The brute 
force method of subdividing the cube into TV 3 sub-cubes and doing the 
sums results in slow convergence because of the Coulomb singularity. 
Using symmetry and algebra the Coulomb singularities can be elim- 
inated. This leads to an accurate numerical algorithm as well as an 
interesting exact result relating the desired interaction energy to three 
other interaction energies, namely those of cubes touching each other 
at a face, at an edge, and at a corner, respectively. As an application 
a simple model illustrating Wigner crystallization is presented. 
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1 Introduction 



There are still many interesting problems involving the electrostatics of cubic 
geometries. These have to do with cubic ionic crystals [1, 2J, with the force 
and potential from cubic charge and mass distributions [3j HJ 16] , and with 
the electric capacitance of the cube [HE]. Here we will discuss the evaluation 
of the electrostatic interaction energy of two coinciding homogeneous cubic 
charge distributions. For unit charge distributions in unit cubes this energy 
is given by, 

c f f dVidV 2 
where, r a = (x a ,y a ,z a ), and 



Vi Jv 2 \r 1 - r 2 \ 



V a = {(x a ,y a ,z a ); 0< x a <l, 0< y a <h 0< z a <l}, a = 1,2. (2) 

This integral arises naturally in the free electron gas theory of conduction 
electrons in metals, see Raimes 0. Its actual value is usually not needed 
for most applications of that theory. In an extension of the theory by Essen 
|10j . however, it is needed. Another application of this type of integral will 
be given below. The value of C can be calculated exactly. Put, 

Mr) = I r-^-r, (3) 
Jv 2 \r - r 2 \ 

for the electrostatic potential energy from a homogeneous cubic charge dis- 
tribution. This potential has been discussed by Waldvogel [3J, by Hummer 
[5], and by Seidov and Skvirsky [6]. Using it we can write (CE]) in the form, 



C= / <i>c{r x )W u (4) 

JVi 

and this makes it possible to find the analytical expression, 

C= - 2 { 2 ^-^-\ | + ln[(V2-l)(2-V3)]}, (5) 

(Seidov and Skvirsky [6J). This evaluates to, 

C w 1.8823126443896601600 (6) 

using twenty digits. Another expression for C in terms of a one dimensional 
integral has been derived by Essen and Nordmark pQ. 

If we displace one of the cubes in (J2]) one unit along the x-axis the integral 
(00) changes into an integral for the interaction energy of cubes with one face 
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touching (see Fig. [T]). Let us call this integral Cf. If we displace one cube 
one unit along both the x and the y-axis we get the integral for cubes with 
an edge in common, call it C e . If we finally displace one of the cubes one 
unit along all three directions of space we get the integral for cubes touching 
at one corner, call it C c . One of the results found below then says that, 

C = Q + C C + l -C c . (7) 

This might be a new result. 

Coulomb interaction energy integrals find one of their main applications 
in Hartree and Hartree-Fock self-consistent field studies of many electron 
systems, see for example Raimes [9]. As an application of the results of this 
paper we use them for crude estimates of the energy of electrons moving in a 
cubic background of smeared out positive charge. In particular we compare 
the energies of delocalized electron states with those of localized states. When 
the density is small the localized states are found to have lower energy. This 
is the phenomenon of Wigner crystallization pj], [12] . 



2 The brute force approach 

Consider two electrons of charge e in a cubic box with edges of length L. 
Assume that both electrons have constant charge density, 

P = e/L 3 , (8) 

in this box. The Coulomb, electrostatic, interaction energy of these charge 
distributions is then: 

4=(£)7 (/ r^-fW ( 9) 

L \L 6 J Jr^V! \Jr 2 GV2 {rx - r 2 \ J 

Here V a , (a = 1, 2) denote the cubic boxes over which the integration vari- 
ables, r a = (x a ,y a ,z a ), take their values. We now introduce units so that 
e = L = 1. The integral can then be expressed in the form, 

^ rx=i ry=i rz=i ru=i rv=i rw=i dx dy dz du dv dw 

=° \j(x-u) 2 + iy-v) 2 + (z-w) 2 ' 



x=0 J y=0 J z=0 J u=0 Jv 



which shows explicitly that this is a six-dimensional integral. 

Nowadays we are spoilt by systems for doing mathematics by computer. 
It is therefore tempting to try these systems whenever some cumbersome 
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integral arises, and frequently they do deliver sensible answers. For the inte- 
gral ffTUj) . however, those that I have tried fail. Brute force can't handle the 
Coulomb singularity. Let us see what happens if we start by dividing each 
cube into N 3 sub-cubes: 

ijk _// \ ^ ~ 1 ^ j — 1 3 ' k — 1 k 1 

(11) 

where the indices, and fc, run from 1 to N. Our integral can then be 
written as the sum, 

TV N 

C= E E Qjfty*, (12) 

ijk=l lmn=l 

over iV 6 terms, integrals over pairs of sub-cubes, 



/ / | V (13) 



For sufficiently large N most integrals are over pairs of spatially separated 
sub-cubes and can be easily approximated. This leads to a brute force ap- 
proach. Fairly large contributions should, however, come from pairs of cubes 
that coincide or touch since they are strongly affected by the singularity. 
Such an approach is clearly clumsy. 



3 Removing the interior singularity 

The awkward singularity occurs only in the interior of those N 3 terms of this 
sum for which the integration sub-cubes are equal. If we thus write, 

N N N ' 

c= E ct Jk + E E cfe> (1 4 ) 

ijk=l ijk=l lmn=\ 

where the terms with all three indices the same (i = l,j = m,k = n) are 
excluded in the double sum, we see that the interior singularities occur in 
the first sum over coinciding sub-cubes. But these integrals are all identical 
and equal to, 

Cs = I (15) 

From formula (114"]) one thus gets, 

JV N ' 

C = N 3 C N + E E (16) 

ijk=l lmn=l 
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Apart from being over a smaller cube, the integral Cjv is essentially like the 
original integral. In fact one easily finds the scaling property, 

C = N rj C N . (17) 

Using this equation (TTBI) becomes, 

(i N N ' 

C =^ + E E (is) 

ijk=l lmn=l 

Solving for C we thus finally have, 

jy2 N N ' 



c 



T E E (19) 



iV 2 

ijk=l lmn=l 



Here the original integral with its singularity has been written as a sum of 
N 6 — N 3 integrals without (interior) singularities. 



4 Approximating the non-diagonal integrals 

The non-singular integrals can be approximated by the product of the two 
cubic volumes divided by the distance between their midpoints. A simple 
calculation gives, 

C lmn = f f dTWa ^ 1 1 m) 

N ' ljk JvWJvlvr ki-r 2 | ~ N s ~ _ ' ; 



yj(i - I) 2 + (J - mf + (k - n)< 



If we introduce the notation, compare equation ffTTj) . 



CgT = iV 5 C^ fc , (21) 

we now have, 

i at at ' 

r~i \ \ r~ilran l r ) r )\ 

°- iV3(AT2_i) 2- 2. °<i* • 

Since the number of terms in the sum grows as iV 6 it is of interest to take 
advantage of symmetries to reduce it as much as possible. Doing this we find 
that, 

19 / at2 N N N 9 N N N 

° " N*(N* _ 1) T" £ + 7V ^. ^ C W + 3 tV ^ h l]k 

\ ' \ Ki m<in<k Ki m<i n<k 



l<i rn<j n<k l<i m<j n<k 



(23) 
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N 




c%-c 


2 


1.899556871 


2 ■ 10~ 2 


4 


1.887296187 


5- 10~ 3 


8 


1.883654361 


1 • lO- 3 


16 


1.882660569 


3 ■ 10" 4 








oo 


1.882312644 






Table 1: This table illustrates the slow convergence of C^. Note that computation 
time goes as N e . 



is an alternative way of writing equation (j22p . Putting, 

AgT = ^(i-Z)2 + (j- m )2 + ( fc -n)2, (24) 

we have that, 

CgT « (25) 
assuming that (z, j, fe) 7^ (l,m,n). If we put this into (122)1 . or (123)) . we get, 

j AT iV ' ^ 

C ~ CiV = /V3C/V2 _ n E E 7\7mn = ( 26 ) 
iV l iV ^ ijk=llmn=l ^ijk 

12 /iV 2 ^ 1 _J_ 2 " * * 1 \ 

iV l^iV Ij y Z ;<i ZA in m<jn<k^jk 6 Ki m<j n<k ^ijk J 

The smaller the box, the smaller the error, so there is hope that this ex- 
pression will converge to the correct value of C when N goes to infinity, i.e. 
that 

C = lim C° . (28) 

The approximation (125)) then immediately gives the following estimate for 
C, when N = 2, 

C « C 2 ° = 1 + + ~ « 1.899556871, (29) 

a value which turns out to be correct to two significant digits. This is en- 
couraging but the convergence for increasing iV is slow, see Table [TJ When 
iV = 16 the error is still 3 ■ 10~ 4 . The exact value is from Eqs. (jSJ) and 
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Figure 1: This figure illustrates the three cases of touching cubes for which 
the integrand is singular on a face, an edge, and a corner, respectively. 



5 Removing the remaining singularity 

We now introduce the symbols, 

W+1,1,1 y-fl,m+l,n+l /~i y"rf+l,m+l,n+l (Qn\ 

U l,l,l — °f' u l,m,n = °cj °i,m,n = u ci \<-> u ) 

for the integrals between adjacent sub-cubes that have a face, an edge, and 
a corner, in common, respectively (see Fig. [T]). These represent the terms in 
the sum ( |23l) that still are affected by the Coulomb singularity. Using this 
notation formula (|23|) for the case N = 2 gives, 

C = Q + C C + ±C C , (31) 
which is the result ([7]) promised in the introduction. 



The integrals (l3"Uj) occur in the sum f[2"3"l) for C the following number of 
times, 



iV f = 6iV 2 (iV- 1), N e = 12N(N — l) 2 , iV c = 8(iV - l) 3 , (32) 
respectively. Let us put, 

Fn = iv 3 (iv 2 - 1) = iv(iv + 1) ' (33) 

N = N 5 (N 2 - 1) ~ N 2 (N + 1)' ^ 

c - ^ 8(iv-i) 2 

and define the two quantities, 

S N = C N — Fn — En~^ — ^ N ~^^ (^) 

and, using this, 

Cat = Sn + F N Cf + E N C C + CnC c . (37) 



iV e = 12N(N 


- D 2 , N c 

/ 5 C 


N f 


6 


' N*(N 2 - 1) 


JV(JV + 1)' 




12(iV- 1) 


W(N 2 - 1) 


7V 2 (iV+ 1)' 


N c 


8(iV- l) 2 


N 3 (N 2 - 1) 


N 3 {N + 1)' 



Clearly Sn is the sum of the terms in (126]) that approximate integrals that do 
not contain singularities (in the interior or on the boundary). We thus have 
that 62 = since for N = 2 all the sub-cubes are in contact. Therefore C N 
is an estimate of the integral C by a sum in which the integrals containing 
surface singularities have been replaced by their (unknown) exact values, 
while the remaining ones are estimated by their inverse distance, Eq. (1251) . 
The function C N obeys both (since <5 2 = 0), 

C\ = C, (38) 

and, 

lim C N = C. (39) 

N— >oo 

Since the 5n are known quantities the assumption that C N = C, in equation 
(j3"7]) . gives for each N an equation in four unknowns (C, Cf, C e , C c ). A system 
of four such equations, 

C - F Nk C { - E Nk C c - C Nk C c = 5 Nk k= 1, 2, 3, 4, (40) 

can thus be solved for these unknowns. Now, each quadruple of numbers 
N\, N2, N3, N4, will give us an estimate of the four integrals. In calculating 
the 5]y k the approximation (1251) has only been used for integrals in which 
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Ni 


iV 2 


iV 3 


N 4 


C 1 


C} 


cl 


cl 


2 


3 


4 


5 


1.882304130 


0.98272866 


0.70632105 


0.57976327 


2 


6 


8 


10 


1.882311519 


0.98306698 


0.70575406 


0.58047142 


2 


11 


15 


19 


1.882312489 


0.98340873 


0.70521257 


0.58107356 


2 


20 


25 


30 


1.882312615 


0.98367876 


0.70479560 


0.58151474 


2 


30 


35 


40 


1.882312641 


0.98390505 


0.70445014 


0.58187235 


2 


44 


50 


56 


1.882312647 


0.98409569 


0.70416088 


0.58216823 



Table 2: The rows of this table illustrate solutions of the system of equations 
(|40|) . C 1 has converged to C in the last rows but the convergence to the other 
three integrals is clearly slow. 



the integrand does not become singular. Obviously one of the numbers 
should always be chosen to be two since then one of the equations of the 
system is exact. 

In Table [2] some results of this approach are shown. After finding four 
different sets of quantities a standard linear equation solver delivers four 
solutions to the linear set of equations. For C this is clearly seen to give 
excellent values. The three other integrals converge much more slowly but 
seem to approach C { « 0.984, C c ps 0.704, and C c ps 0.582, respectively. 



6 Electrons in a homogeneous cube 

Here we will use crude estimates of the Hartree energy [9] of electrons that 
move in a cube of homogeneous positive charge density. Using this crude 
theory we will investigate whether the electrons tend to delocalize in the 
cube or if a state with localized electrons has lower energy. 

We assume that the electrons either are delocalized in the cube and have 
constant charge density in the cube or that they localize in one octant of 
the cube and have constant charge density there. This means that we can 
treat either 8 electrons or 8 electron pairs. With these assumptions the 
electrostatic interaction energy can be found from the results above. The 
kinetic energy is estimated essentially by means of the uncertainty principle 
and the Pauli exclusion principle. 

We start with 8 electrons in delocalized states. They are assumed to move 
in a cube of side L and positive charge 8e. The energy is then the sum of 
the kinetic energy, 



2 

/s " = 2^ 



3 2 1 

2— + 6 — + 



L 2 \L 2 {L/ 2) 



2 



(41) 
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and the electrostatic energy, 



In the kinetic energy the two first electrons are assumed delocalized over the 
cube without nodes in the wave function. The remaining six must then go 
into the three degenerate states with one node. The first term in the elec- 
trostatic energy is the self energy of the positive background. Then follows 
the attraction between the background and the eight delocalized electrons. 
The final term is the sum of the 8-7/2 electron-electron pair repulsion terms. 
Simplifying this gives the total energy 

ft 2 42 „ „ 4 

E ™ = ^v- Ce T < 43 > 

The energy of this closed shell delocalized state should now be compared 
to the energy of the ferromagnetic localized state with the electrons in one 
corner each. We find, 

T - k \ 3 - h * % (44) 

for the kinetic energy since all 8 electrons now sit in cubes (octants) of side 
L/2. They are however alone in their corners (octants) so the Pauli principle 
is automatically obeyed. The electrostatic energy becomes 

Vsl = l^L~° ~ *T/2 {C + 3Q + 3Cc + Cc) + lT (3C<f + 3Cc + ° c) ■ (45) 
Simplification of this using Eq. (131 j) gives the total energy 

E - = ^- Ce L- (46) 

If we introduce atomic units (h — e — m — 1) so that length is measured in 
units of the Bohr radius we can plot the two energy curves, 

21 4 

E 8d = Jji-C^ (47) 
Esi = §~Cj, (48) 

and get the results of Fig. [2j 
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Figure 2: The energies of eight electrons in a cube of side L as function 
of L for the localized and delocalized cases, respectively. For L > 3.59 the 
localized state has lower energy. 



Finally we give the corresponding results for 16 electrons sharing orbitals 
pairwise. In the electrostatic energy one can then essentially change the 
particle charge e to 2e and add the contributions from the repulsion within 
the pairs. This gives the two curves, 

Em = %-c\, (49) 
E m = ™- £ (50) 

for the delocalized and localized energies respectively. These curves are plot- 
ted in Fig. El 

One notes that the localized states always have lower electrostatic energy 
simply because in these states the electrons are better at avoiding each other. 
For small L- values the delocalized states always have lower energy because of 
the uncertainty principle. The curves in these plots resemble those of Wigner 
[TT1 [T2] who predicted that localization gives lower energy in metals at low 
densities. This phenomenon is called Wigner crystallization. 



7 Conclusions 

I am not aware of any comparable study of the electrostatic interaction en- 
ergies of homogeneous cubic charge distributions. The algebraic and com- 
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Figure 3: The energies of eight electron pairs in a cube of side L as function 
of L for the localized and delocalized cases, respectively. For L > 2.49 the 
localized state has lower energy. 

binatoric tricks used to eliminate the Coulomb singularities in the integrals 
seem partly new, as well as the result of Eq. (J7j). It is possible that these 
ideas can be generalized to more general integration problems involving the 
Coulomb singularity. It is a further bonus that these insights into the elec- 
trostatics of cubes and their sub-cubes can be used to make simple estimates 
for the Hartree energy of electrons distributed in cubes in different ways. 
Such simple model systems are of value for the qualitative understanding of 
more complex systems. 
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